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We describe how density-matrix renormalization group (DMRG) can be used to solve 
the full-CI problem in quantum chemistry. As an illustration of the potential of this 
method, we apply it to a paramagnetic molecule. In particular, we show the effect of 
various basis set, the scaling as the fourth power of the size of the problem, and compare 
the DMRG with other methods. 



I. INTRODUCTION 

First-principle quantum chemistry is employed successfully to obtain thermochemical data; molecular 
structures; force fields and frequencies; assignments of NMR-, photoelectron-, E.S.R-, and UV-spectra; 
transition state structures as well as activation barriers; dipole moments and other one- or two-electron 
properties. Two routes of calculation are available: (i) ab initio Hartree-Fock (HF) and post-HF meth- 
ods; and (ii) density functional theory (DFT). Though both approaches are rigorous, the former one 
necessitates lengthy configuration interaction (CI) treatment to account for electron correlation; whereas 
the latter one crucially depends upon the quest for accurate exchange and correlation functionals. The 
recently acquired popularity of DFT stems in large measure from its computational efficiency, allowing it 
to treat medium to large size molecules at a fraction of the time required for HF or post-HF calculations. 
More importantly, expectation values derived from DFT are, in most cases, better in line with experiment 
than results obtained from HF calculations. This is particularly the case for systems involving transition 
metals. Nevertheless, if one wants to achieve experimental accuracy for small polyatomic molecules, the 
method reaches its limits. On the other hand, post-HF overcomes these limits and goes even beyond 
experimental accuracy, e.g. in case of H2. The drawback is of course its very high cost in computational 
power. However these very accurate calculations of small systems may provide the best route to obtaining 
more accurate exchange-correlation potentials using the constrained search algorithm of Levy ]IJ . 

Keeping the discussion to ab initio (HF and post-HF) quantum chemistry, let us now consider the state 
of the art in this topic. We now know how to do very large calculations, using the direct methodology 
0. We can also manage to work with good basis sets for such calculations, although it is considered 
that 6-31G* are not good enough, and probably something nearer to TZ2P is required for definitive 
SCF calculations. Beyond SCF there are major difficulties, all associated with trying to more accurately 
represent the electron-electron cusp. We know from the work of Kutzelnigg |3j, that the convergence 
of this problem is very slow, something like [I + ^) -4 , with I the orbital angular momentum quantum 
number. This means that very large basis sets are required for correlated calculations. We know that 
it is more important to include d and / basis functions than to improve the methodology (e.g. 6-31G* 
basis are not appropriate for correlated studies). We also know that the raw cost of the main correlated 
methods, MP2, MP3, MP4, CISD, QCISD, and QCISD(T) increase with powers of the size of the problem 
as 5, 6, 7, 6, 6, 7, respectively. We recognize that QCISD(T), with large basis sets, can make predictions 
for small molecules which are of "chemical accuracy" (e.g. the G2 quantum chemistry of J. Pople at al.). 
Considering "coupled cluster" technology, although this method is very promising, we have to realize 
that it is not yet possible to contemplate calculations with 1000 basis functions. The rapid progress that 
computer companies have made in the development of their hardware, such as more memory, vectorization 
compilers, and parallel machines will not cure these basis set problems. 

Recently, one of us (S.W.) and R. L. Martin have shown that with the density matrix renormalization 
group (DMRG), one can calculate "exactly" (full CI) the electronic structure of N interacting electrons 
in a field of M nuclei Q. By "exactly" we mean that the correctness of the treatment of correlations 
increases very rapidly and systematically with the calculational effort, so that one is essentially only 



limited by the accuracy of representing the Hamiltonian in a finite basis set. Moreover, these authors 
showed that the method theoretically scales as 0(N ) and hence the calculation of much larger molecules 
will be feasible. More precisely, this ./V 4 dependence originates from the summation over four indices in 
the second quantized form of our Hamiltonian (Eq. 13). Hence, it will be possible to exploit the short 
sightedness of interatomic interactions in the future e.g. in using localized basis function (cf. vide infra). 
Thus we expect to achieve linear scaling, O(N), for large linear molecules when teraflop CPU's and peta 
bytes storage media will become routinely available to us. Moreover, we would like to point out that the 
method at hand is well suited for CAS calculations. Indeed, since the basis functions used will usually 
be Hartree-Fock or Kohn-Sham eigenfunctions, we can easily use the corresponding eigenvalues to setup 
an active space by discarding inner shell orbitals (in keeping their occupation frozen) or in eliminating 
those with very large energies. 

This paper is organized as follows. In Sec. II, we recall the physical principles of the DMRG. Then we 
reformulated, in a more chemical language, the basic algorithm for molecules in Sec. III. Finally, in Sec. 
IV we apply the technique on a paramagnetic molecule and discuss the performance of the method as a 
function of the basis set used. 



II. PRINCIPLES OF THE DENSITY MATRIX RENORMALIZATION GROUP (DMRG) 

A. The Renormalization Group 

The idea of the Renormalization Group (RG) is to apply a transformation to the Hamiltonian which 
eliminates degrees of freedom that are unimportant for the description of the states in the energy range 
we are interested in. For example, if we are interested in the low energy states of a system with a energy 
cutoff A, one integrates out modes with energy 

A-SA<uj <A (1) 

where 5A is a small energy interval. Then we rescale the parameters of the new system so that it 
reproduces the previous one. In other words, given a Hamiltonian H of a system with N variables, a 
renormalization transformation is a mapping in the Hamiltonian space which maps H to H' , 

H' = R b {H). (2) 

This new Hamiltonian describes a system with N' = N/b < N variables. The RG transformation must 
be unitary, i.e. it has to preserve the partition function (Z = Tr e~ H / kT ), 

Z N , [H'\ = Z N [H] . (3) 

This set of transformations is called a group because the transformations must be associative: 

R h , (Rb(H)) = R b , b (H). (4) 

An exact transformation is usually not possible. For example, Wilson and Fisher worked in Fourier 
space and used a perturbative scheme in order to analytically solve this problem j^] . Another approach 
is to consider the Hamiltonian on a lattice and to work in real space. A RG transformation could be 
the elimination of every second site and then the redefinition of the Hamiltonian parameters in order to 
reproduce the original problem. 

A numerical implementation for a quantum mechanical problem was given by Wilson who succeeded in 
constructing a non-perturbative RG transformation for the Kondo model ||. The same transformation 
has since then been applied to a wide range of similar quantum impurity models with equal success 
0. On the other hand, it has proven to be more difficult to develop similar RG transformation for 
quantum lattice models. These difficulties with RG approaches in real space provided the motivation for 
the development of DMRG. 



2 



B. Wilson's RG 



Let H be a one-dimensional Hamiltonian describing an interacting system of electrons on a lattice with 
L sites labeled i. Each sites has one (Wannier) orbital which can be in 4 different states : |0), | t)> | 4) 
and | ][). The dimension of the Hilbert space for N-\ up-spin electrons and JV| down-spin electrons is 
thus 



dimH = 




Which is, for L — 100 and N* = N\ = 50, dim7i = 10 58 . This is of course intractable numerically. The 
idea is to obtain the low energy eigenstates of that system by keeping only a small number of states, e.g. 
1000. 

For the sake of clarity we will first directly apply Wilson's scheme to quantum lattice systems and then 
show why it does not work with a simple example. We will describe next a RG step. 
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FIG. 1. Wilson's blocking scheme for real-space renormalization group. 



Let Bi be a block describing the first I sites for which we only keep m states to describe the Hamiltonian 
of that part of the system. The internal Hamiltonian Hg is then represented by an m x m matrix. One 
adds to this block (see Fig. another block By describing £' sites with internal Hamiltonian Hy . The 
Hamiltonian of the new block Bg + y has dimension mm! and can be written as 

H l+ y = H t ® 1 + 1® Hy +^c a A a ®B a . (6) 

Here, the last term represents the interactions between Bg and By, where the A a are operators in Bi 
and the B a are in By (See Section III.C.) The RG idea is then to reduce the size of this many body basis 
from mm! to m. For this purpose, in Wilson's original formulation, one diagonalizes Hg + y and keeps 
only the m states with the lowest energy eigenvalues, since we seek the low energy states. 

This method has been directly applied to quantum lattice systems || but the spectra become inaccurate 
after a few iterations. The problem originates from the treatment of the boundary conditions. To illustrate 
this, we will consider the problem of a free particle on a lattice of size L with open boundary conditions 
or, in other words, the discretized version of a free particle in a box. The eigenstates of this system are 

Xn(x) ~ sin (Jy x ) ( 7 ) 

where n is a positive integer. We now consider the addition of two blocks of size L. In the RG procedure, 
the lowest few eigenstates of system of length L are combined to form the low-lying eigenstates of the 
system of length 2L. The lowest eigenstates of two systems of length L and of a single system of length 
2L are shown in Fig. |[ One sees clearly that a combination of the ground states of the two small systems 
(size L) cannot reproduce the ground state of the large system (size 2L). The kink in the superposition 
of the two ground states of the small systems can only be removed if one keeps almost all the eigenstates 
of the smaller systems. 
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FIG. 2. The lowest eigenstates of two 16-sites system (open circles) and 32-sites system (full squares) for a 
one- dimensional free particle with open boundary conditions. 

This example shows that the treatment of the boundaries of the block is crucial. White and Noack B 
formulated two types of RG procedures which solve this problem. Both are based on choosing a new basis 
for H 2 l which is not the set of eigenstates of Hl . In their first procedure, the combination of boundary 
conditions method, the new basis is obtained from the low-lying eigenstates of several different block with 
different boundary conditions at the edge of the block. In the second type of procedure they proposed, 
the superblock method, the new basis for H^l exploits the idea that it will eventually be used to make up 
parts of a larger system. These two methods were applied to the free particle in a box and both yielded 
accurate results. While the first method is not so well suited for many particle problems, the second one 
necessitates an optimal projection from the superblock. This gave rise to the DMRG scheme which is 
described in the next section. 



C. Density-Matrix approach 

The fundamental difficulty in the standard approach lies in the choice of the lowest energy eigenstates 
of Hi + i> as the new truncated basis. Since Hn+i* contains no connections to the rest of the lattice, its 
eigenstates exhibit inappropriate features at the block boundaries. One way of avoiding this problem is 
to diagonalize a larger system (a "superblock") which contains the new block Bi + gi. The idea is that 
fluctuations in the rest of the superblock (the "environment" ) lead to a more appropriate treatment of 
the boundary of Bg + ii. The wave function for the superblock is projected onto Bi + ii producing a set of 
states of Bz + ?j which are then retained. 

Let {\i)} be a complete set of states of Bt+ti and {|j}} be the states of the rest of the superblock. We 
can then write 

W) = Y.^Mf)- (8) 

In Ref. |l0| one of us (SW) shows that the optimal states to be kept are the eigenvectors of the reduced 
density matrix of the superblock, 

Pij = ^iptkipjk- (9) 

k 

Let us denote by u v the eigenvalues of p in decreasing order. The m optimal states we seek, u v , are the 
eigenstates of p with the largest eigenvalues. Each u> v represents the probability of the system being in 
the state u v ', hence J2u UJl/ ~ ^ ne deviation of 

m 

i/=i 

from unity measures the accuracy of the truncation of the m states (typically, 1 — P m < 10 -6 ). 
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1. Density-Matrix algorithm 



In practice we take one block to be just a single site. Fig. g shows the superblock configuration that 
we use. We adopt the notation Be • •Bgi for this configuration, where Be represents a block composed of 
i sites, Bp is a block of length • represents a single site, and the length of the superblock is therefore 
L = £ + £' + 2. Hence, the new block Be+i is formed by the left block plus a single site. 




FIG. 3. The configuration of blocks used for the density matrix calculations. 



In general, DMRG algorithms can be divided into two classes according to how the environment block 
Be/ is chosen. In the infinite system algorithm, B^ is chosen to be Bf, the reflection of Be- Therefore 
the superblock grows as Bg is built up. For the finite system algorithm, Be' is usually chosen so that the 
length of the superblock is L. It begins with the use of the infinite system algorithm (see Fig. |j) for 
■§ — 1 steps, such that the final superblock obtained is of size L. 



□ ooD 
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L 

FIG. 4. Schematic infinite system algorithm. 

That is, we start with a four site chain B\ • »Bi . Then we calculate the density matrix and construct 
an effective Hamiltonian for Bi- In the second step, we diagonalize Bi • »B^ where we have formed B^ 
by reflecting B^- We continue in this manner, diagonalizing the configuration Be • »Bf , setting B' e , 1 
equal to Bi», and then using Be+i and its reflection in the next step. Once the system Bl_ 1 • »5f_ 

2 2 ~ 1 

is used to form Bl, the next step is to form Bl +1 . This system and all the other superblocks to follow 
will involve exactly L sites. We continue to form the other blocks up to size L — 3 using the superblock 
Be • •-B|*_^_2 to build Be+i- This sequence of steps corresponds to the first iteration of the finite system 
algorithm. 

noo i i 

I |o o| I 



L 

FIG. 5. Schematic finite system algorithm 
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The second and subsequent iterations use superblocks of fixed length L (see Fig. g). We start by 
diagonalizing B\ • »B^_ 3 forming Bi up to -Bl-3 • where we use the blocks obtained from the 
previous iteration to build the reflected blocks Bf . On the very last iteration we stop after diagonalizing 
Bl_ 1 • •■B? _, and then use the wave function of the L sites system to measure ground-state properties, 

2 2 1 

e.g. correlation functions. At each sweep (left-to-right and right-to-left iteration) we increase the number 
of states m used to represent the Hamiltonian and other operators inside the block. 

Since we truncate the basis, the algorithm is variational and in practice the calculated value for the 
ground-state energy obeys the empirical law 

E(m) = E + Ae~ m / m * (11) 

where Eq is the targeted energy; A and m* are adjustable parameters. 

An important improvement in efficiency can be made by keeping track of the wave function in the 
process of adding a site [jfl|. We can then construct a better initial wave function for the Davidson 
diagonalization of the superblock, thus reducing the number of steps needed to obtain convergence. 



III. QUANTUM CHEMISTRY WITH THE DMRG 

A. The non-relativistic time-independent Schrodinger equation 

We will now show how the DMRG can be applied to ab initio quantum chemistry . Essentially the 
DMRG will yield both the energy and the wave function of the ground state and of some excited states. 
More precisely, we consider a molecule with N interacting electrons moving in the field of M fixed nuclei. 
The Hamiltonian describing this system is 

N N M 2 2 

i=l i=l fc=l 11 K| i,j ' Jl 

and we will use atomic units throughout this paper. 



B. Second quantized form of the Hamiltonian 



To perform calculations, this Hamiltonian is represented in a finite basis of K functions. One chooses 
usually atomic orbitals Xi with i = l,...,K which are described, for example, by the product of a few 
Gaussians and spherical harmonics. This set expands a smaller (finite size) Hilbert space and is not 
orthogonal. The first step is to orthogonalize these functions to be able to represent the Hamiltonian 
( |l2| ) in second quantization. In practice one can perform a Hartree-Fock calculation in order to generate 
a set of basis functions. In doing so, we have a good (low energy) starting point for the calculation of the 
correlation energy. 

Consider a set of K molecular orbitals (MO), {<fi}, i — 1,...,K, which we reorder with increasing 
energy. We mainly used Hartree-Fock molecular orbitals, but we do not have to restrict ourselves to 



them, and we have also used Kohn-Sham orbitals (see discussion in IV B). In case of closed shell systems, 
an approximate ground state would be when the first N/2 orbitals are occupied by 2 electrons and the 
remainders (virtual orbitals) are empty. These MO can be obtained from the atomic orbitals by the 
well-known Linear Combination of Atomic Orbitals - Self Consistent Field (LCAO-SCF) scheme jjjj. 
One uses then the second quantization formalism to reexpress the Hamiltonian ([12]). An introduction 
to that technique can be found, e.g., in Ref. p3| . The creation (annihilation) operator c\ a (c icr ) creates 
(annihilates) an electron with spin a in the MO </?j. The Hamiltonian we shall use from now on reads 
then 

H = YJ2 H i3 c L c i° + ^Yj Gi W c \A°' c *°' c te- ( 13 ) 

i.j (7 i,j,k,l a,(7' 

The sum is over all orbitals and spin a. The integrals 
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Ha 



dr(p*(r) 



1 M 7 

9 ^ r - 



fc=l 



^e 2 
r - r k | 



(14) 



represent the elements of the core Hamiltonian (kinetic energy and nuclear attraction) , while 

G m = J J drdrV*(r)^(r')^-^^.(r') W (r) 

are the interelectronic repulsion integrals. Using the occupation number formalism, a Slater determinant 
ijj of N electrons can be written as 



(15) 



^ C i 1 ai C i 2 (J2 



(16) 



where |0) is the vacuum state. As an example, the singlet Hartree-Fock eigenstate for a closed shell 
system is 



N/2^N/2l 



|0). 



(17) 



A full-CI calculation is then the exact diagonalization of (113) in this occupation number representation 
of the Hilbert space basis. Since the size of the matrix of ( |13| ) is 0(N\), e.g. 10 spin-up electrons and 10 
spin-down electrons with 40 orbitals has a dimension of approximatively 10 18 , one cannot use the Lanczos 
|l4| or Davidson (ll| algorithm to solve this problem. Instead we use the DMRG to obtain the ground 
state and some excited states of the Hamiltonian ( fl3| ) for small to medium sized polyatomic molecules. 

Since the DMRG was initially developed for one-dimensional quantum lattices, a molecular wave func- 
tion (such as the Hartree-Fock MO) becomes a site which can be empty, singly occupied with an up or 
down-spin electron or doubly occupied. Hence a molecule described by K orbitals becomes a lattice with 
L = K sites and the Hamiltonian given in Eq. (|l3|) . As input for the DMRG we need the matrix elements 
Hij and Gy^. Since our Hamiltonian conserves the total number of particles and the total spin, we also 
specify the number of electrons N and the projection of the total spin on a quantization axis S z in order 
to reduce the dimension of the Hilbert space. 



C. One DMRG step 

We want to describe here explicitly the major part of the DMRG applied to quantum chemistry, a 
renormalization step with the density matrix. We use the finite-size algorithm of the DMRG and target 
only the ground state. In a warmup procedure, we construct an approximated ground state (e.g. Hartree- 
Fock). Then we start the algorithm following the description in Ref. (lj|. Let us consider in detail the 
addition of an arbitrary orbital (site) . We already have built a subspace Sg which describes partially the 
electronic structure of the molecule using the first £ MOs. In the original DMRG formalism this subspace 
spanned by £ basis functions, which were Wannier functions located on each site, was called a block. To 
avoid possible confusion with quantum lattice concepts, we prefer here to use the word "subspace" . For 
the description of that subspace we keep only m states. These states are used to represent i) the part of 
the ground state in that subspace; ii) the internal part of the Hamiltonian, i.e. terms that only involve 
indices which belong to the subspace (i,j,k,l < £); and iii) terms that will be used in constructing the 

interactions between the subspaces, such as c\ a , [ci a ] , cj^cj^ 

We add then the (£ + l) th orbital (site) to the current subspace (see Fig. ||) and construct the internal 
Hamiltonian for the new subspace Se+i and for all other operators needed as well. 



• 









l+l 



FIG. 6. Schematic view of the left subspace with one orbital added. 
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The Hamiltonian reads 



H t+ i = H e <g> 1 + 1 (g) F x + 



[Q+1,<t]i 



(18) 
(19) 



l,j<£+l cr,a' 



where [A]i denotes the matrix of the operator A represented in the basis of the to most significant states 
of the subspace Si and [A] i denotes the matrix of the operator A represented in the Hilbert space of one 
orbital. For example, if the order of the basis set is {|0), | |), | J.}, | Tl}}> where | fl) is defined as cjcj|0), 
the matrix for the creation operator of an up-spin electron for the MO I + 1 is 



,t 





10 



10 



(20) 



The new subspace is now described by 4m states. We will illustrate now how to reduce this number to 
to. To this purpose we need to construct the environment which involves all the residual orbitals. In 
case of a molecule, we can not reflect the current subspace (block) as in case of a quantum lattice system 
since there is, in general, no equivalent reflection symmetry. We apply the same procedure for the right 
subspace Sl-i-2- In Fig. 0we show the construction of the left and right subspaces. 



>L-l-2 



FIG. 7. Schematic view of both subspaces. 

At the end we will have two subspaces Sf+i and Sl-i-i] the first describing partially the electronic 
structure of the molecule within the first I + 1 orbitals and the second one, the rest of the molecule with 
the remaining L — £ — 1 orbitals. In order to retain only the to states of Si+i which best represents 
the whole molecule, we construct the full molecule by adding the two subspaces together. We build a 
representation of the full Hilbert space by adding the two subspaces together, which means taking the 
direct product of the two subspaces. We then compute the ground state of the whole molecule using 
the Davidson algorithm. |15| We then construct the projected density matrix. The full Hamiltonian 
could be constructed in the same way as when adding two blocks together, Eq. (19), but to be more 
efficient, we use the following procedure. In the Davidson algorithm, the time consuming part is the 
H\ty) multiplication. The matrix of the molecular Hamiltonian can be written in the general form as 



H 



molecule 



ij-i'j 1 



A%B% 



(21) 



where the index i and i' stands for states in the left subspace and for the right subspace. The sum 
on a is over all terms of the Hamiltonian. Then the product jy m °l ecu l e | \jr ^ can be written as 



i',3 1 



H 



molecule 



(22) 



The last sum is performed first as a matrix multiplication of B a with ^ T , to form a temporary matrix 
C%. Then a matrix multiplication of A a by [C a ] T yields a partial term, which is added into the resulting 
vector, giving a sum on a. Using the Davidson algorithm one obtains the ground state and some 
excited states, if desired. 
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Then we construct the reduced density matrix tracing out the environment components 

3 

One diagonalizcs p, obtaining 4to eigenvectors and the corresponding eigenvalues uj v . Each uj v repre- 
sents the probability of the system being in the state <p v , with 

4m 

E w " = 1 - ( 24 ) 

v=l 

One chooses the to states with the highest w„ value. The discarded weight 

m 

l-$^w„ (25) 
i/=i 

measures the accuracy of the truncation to to states. 

The reader may notice a similarity between this procedure and the construction of natural orbitals. 
However, a crucial difference is that the density matrix here is a many particle density matrix, with an 
exponentially large dimension, not a single particle density matrix. The resulting basis is a many particle 
basis, with each state representing a complicated sum over large numbers of configurations. 

When the new basis is formed (to elements), the operators must be updated. The eigenstates <j> v can 
be written in the form 



4>ij (26) 

where i denotes a state in the subspace St and state j corresponding to the (£ + l) th orbital. We store 
them in a to x 4m matrix O. Thus each operator A that is needed is replaced by 

A = OAO T . (27) 

Adding the (£ + l) th orbital to the subspace Sg corresponds to one DMRG step. There are a total 
(K — 2) DMRG steps to zip through all K orbitals, and another (K — 2) DMRG steps to zip back to 
the beginning, with the roles of the left and right subspaces switched. This is a complete DMRG sweep, 
keeping always m states per RG step. After one iteration we increase the number of significant states 
to, e.g. by doubling it, and start again using the preceding result for Sl-(-2 to describe environment 
(i.e. the complementary space). Since we truncate the basis of the full Hilbert space but calculate the 
Hamiltonian exactly within this basis, this algorithm is variational. In the limit where we keep enough 
states to to solve the full-CI problem, this algorithm is also clearly size consistent. 

At the end, we obtain the energy and a representation of the ground state from which we can compute 
the density, the natural orbitals, the two-electron density matrix, etc. 



D. Illustration 

In order to illustrate the DMRG procedure, let us consider a simple example, the methane molecule 
with Td symmetry (CH4). We follow the procedure outlined above : 

1. Choose a minimal basis : STO-3G (9 orbitals). 

2. Perform a Hartree-Fock calculation obtaining 9 MO and the Hartree-Fock electronic energy E = 
-53.248968... 

3. Represent the Hamiltonian in this basis, i.e. calculate all Hij and Gijki- 

4. Apply the DMRG to obtain the ground state of the molecule (sec Appendix for some technical 
points). 



9 



In Fig. § we show the ground-state energy E as a function of the number of significant states m retained 
in each DMRG iteration. We also show results obtained by the CAS method. In CAS(v,[i) we consider 
the complete active space configuration interaction of fj, electrons for v orbitals. We see that DMRG is 
variational and converges to the CAS(10,9) result which is full CI. 



-53.21 
-53.23 
-53.25 • 
-53.27 
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-53.33 
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CAS(10,9) 




#DMRG 
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FIG. 8. Ground-state energy E of CH4 as a function of the number of states retained m at each DMRG 
iteration. The lines are the Hartree-Fock and CAS results. 



In Table | we also show the occupation numbers 

(rii) = (4 T c lT ) + (c^Cii) 



(28) 



of the MO as a function of the progress of calculation. For instance it can be seen that the occupation 
of symmetry equivalent orbitals (t 2i ) becomes equal at the last iteration. DMRG achieves the accuracy 
of CAS(8,8) (occupation of la\ orbital frozen) when m > 32. The analysis of the occupation numbers 
shows the importance of the higher energy states in the description of the ground state of this typical 
closed-shell system. A significant occupation of the orbitals usually unoccupied in the HF procedure is 
observed. 



Irreducible 


Hartrec Fock 


(rii) after 


(rij) after 




{rii) after 


representation 


eigenvalue 


warmup 


first iteration 




last iteration 


lai 


-11.02969 


2.000000 


1.999997 




1.999965 


2ai 


-0.91210 


2.000000 


1.998077 




1.985416 


u 2l 


-0.52050 


2.000000 


1.996863 




1.976450 


u 22 


-0.52050 


2.000000 


1.992911 




1.976450 




-0.52050 


2.000000 


1.988980 




1.976450 


2t 2l 


0.71922 


0.000000 


0.021086 




0.021952 


2t 22 


0.71922 


0.000000 


0.005932 




0.021952 


2h 3 


0.71922 


0.000000 


0.010275 




0.021951 


3ai 


0.76106 


0.000000 


0.009580 




0.019412 



TABLE I. Occupation numbers (n,:) of the MO as a function of the progress of the calculation. 
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IV. RESULTS AND DISCUSSION 



In order to demonstrate the performance of the DMRG when applied to quantum chemical problems, 
we have chosen a simple magnetic system which has already been the subject of extensive post-HF and 
DFT studies. In fact, computing the energy difference between multiplets of different spins is a difficult 
task for standard computational approaches. This is due to the very high accuracy required for the 
calculation of the ground- and excited-state energies and the intrinsic need of taking into account the 
mixing with states of higher energy (configuration interaction). For this reason, we believe that this 
calculation can be a useful benchmark to test the performance of DMRG in quantum chemistry. The 
crucial features of any new quantum chemical method are scaling and convergence. In the following we 
address this point, both when starting from localized and non- localized Hartree-Fock MO as well as from 
Kohn-Sham orbitals. 



A. HHeH 

The linear HHeH molecule (Fig. ^) was chosen as a typical example for a model magnetic system. In 
this simple case the two paramagnetic centers, carrying each an unpaired spin, are the hydrogen atoms. 
They are linked via a diamagnetic bridge constituted by the He atom. The spins of the two paramagnetic 
electrons can be parallel or antiparallel yielding two different spin states, namely a singlet (S = 0) and a 
triplet (S = 1). The difference in energy between these two spin states as a function of the H-He distance, 
has been subject of both post-HF and DFT calculations [ p7|Jl8| , and shows the importance of electron 
correlation. 




FIG. 9. HHeH: A schematic sketch 

In Table II the results found in literature [|l7| are compared to the exact results obtained with both 
full-CI approach and DMRG using the same basis set (6-311G** ^lf). This basis set contains a total of 
18 orbitals. The following notation is used: CAS(^,/i) = complete active space configuration interaction 
considering v electrons and fi orbitals; QCISD(T) = quadratic CI including single and double excita- 
tions and also triple excitations perturbatively; BS = broken symmetry p0| ; SD = single determinant 
p9| ; LDA = local density approximation for exchange correlation functional in the Vosko-Wilk-Nusair 
parameterization [^8| ; BP = generalized gradient approximation for exchange and correlation functional 
using the Becke parameterization |^9| . The data obtained with post-HF methods show that a large part 
of electron correlation is already included at a modest CAS level. Both DF based techniques, the single 
determinant method |h| (SD) and the broken symmetry approach 
terms of computational resources, are just adequate to yield qualitative results. 



Method 


Distance=1.25 A 


CAS(2,2) 


4204 


CAS(4,3) 


4294 


QCISD(T) 


4298 


BS-LDA 


12432 


BS-BP 


10529 


SD-LDA 


9050 


SD-BP 


7799 


Full-CI 


4860 


DMRG 


4859 



TABLE II. HHeH Singlet-triplet energy gap (in cm ') as a function of the H-He distance for various methods. 
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In Fig. [h] the singlet-triplet energy gap as a function of m is reported for a H-He distance of 1.25 A. 
For m > 64 the DMRG procedure has converged to the exact result (full CI). Results superior to those 
obtained with DFT are observed already when m w 32. The same level of accuracy as with CAS(4,3) is 
also reached for this m value. 
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FIG. 10. Singlet-triplet energy gap of HHeH as a function of m for a H-He distance of 1.25 A 



B. Choice of the basis 



The results obtained with DMRG when applied to quantum lattice systems [|22j suggests that the 
use of localized basis sets, which reduce the interactions between far distant neighbors, will enhance 
convergence. Therefore we attempted to localize the starting MO basis set using standard localization 
techniques. Amongst the classical quantum chemical localization techniques ]2^-|25[|, we have chosen to 
apply the procedure of Mezey and Pipek |26| due to the fact that just overlap integrals are needed to 
perform the localization. This procedure is based on the minimization of a mean localization measure of 
the occupied molecular orbitals D defined as 



N 



D-^N-^d- 1 



i=l 



with 



(29) 



(30) 



fe=i 



and where k are the different nuclei of the molecule, i the molecular orbitals, and Q\ are the gross atomic 
Mulliken population of the orbital i. Rather than localizing all the MO, we only localized the subset of 
the virtual MO (not occupied at a HF level). The two subsets, formed by the occupied non-localized MO 
and the virtual localized MO, are still orthonormal. In fact, the localization of the whole set of MO would 
yield a substantial increase of the starting energy of the DMRG procedure, while localizing the virtual 
orbitals only does not. Recalling the empirical dependence of the electronic energy upon m (Eq. ^), we 
show in Fig. [H], E(m) — Eo versus m to illustrate the convergence of the DRMG procedure using localized 
and non- localized Hartree-Fock MO. Unfortunately, we observe that the localization procedure does not 
yield a faster convergence. In fact, it should be stressed that this system is not the most adequate to 
test the impact of localization. Indeed, there are only a few electrons involved and the spatial extent of 
the molecule is rather small. Nevertheless, we expect that this localization procedure should improve the 
convergence of the DMRG for larger systems, e.g. linear chain like molecules for which most of the Vijki 
terms will vanish or become negligible when using localized basis functions. 
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FIG. 11. Rescaled energy E — Eo as a function of m using localized and non localized HFMO on a lin-log scale. 
The dashed lines are just guides for the eyes. 

Another feature we investigated was the possibility to use Kohn-Sham orbitals as basis functions. We 
performed DMRG calculation starting from Hartree-Fock orbitals as well as from Kohn-Sham orbitals 
obtained using both LDA [^8| and GGA |27|| fu nctional. The results obtained for the singlet state with 
a H-He distance of 1.25 A are shown in Fig. |12|. We see that all three basis functions have essentially the 
same convergence rate. 




100 200 

m 



FIG. 12. Rescaled energy E — Eo as a function of m using HF, LDA and GGA orbitals, on a lin-log scale. 

C. Scaling of the method 

In order to test the numerical effort spent with DMRG we performed several tests on the HHeH singlet 
state at a distance of 1.25 A, either using different basis sets (varying K at constant m) or keeping the 
basis set (varying m at constant K). The aim of these calculations was twofold. First, to test the scaling 
of DMRG H and second, to demonstrate its superiority with respect to standard post-HF approaches 
(full CI, QCISD or QCISD(T)). The main disadvantage of post-HF solutions (apart from their scaling) 
is the lack of at least one of the main feature of DMRG: size consistency and variational character. The 
only theoretical approach which fulfills all these requirements is presently full CI. For example, CID or 
QCISD methods are variational but not size consistent. On the other hand, all perturbative techniques 
(MPn), whose scaling are comparable with the one of DMRG, are not variational and practically feasible 
only up to n = 4 (inclusion of quadruple excitations). 
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The theoretical scaling of DMRG (neglecting the cost of the preliminary Hartree-Fock procedure) is 
0{K A m 2 ). In our case K is relatively small (K max = 30), while m can be an order of magnitude bigger 
than K (m ma _ x = 264). In Fig. [l3] the plots of the CPU time versus K at fixed m and versus m at 
fixed K are shown. The slope, obtained by linear regression, observed in part a) is in agreement with the 
theoretical predictions, i.e. 0(K A ). While the slope in part b) is lower that predicted since we use highly 
optimized BLAS to perform the matrix operations (cf Appendix). Thus the scaling of DMRG appears 
to be competitive even with QCISD, which scales as 0(K 7 ), while yielding much better results. 
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FIG. 13. a) CPU time as a function of A" on a log-log scale; b) CPU time as a function of m on a log-log 
scale. Both plots are for HHeH. 
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APPENDIX A: IMPLEMENTATION OF THE METHOD 



The program we use is written in C++, an object-oriented language. We have to keep track of all sorts 
of matrices connecting different subspaces with different quantum numbers, which is conveniently done 
using classes and overloading operators. In addition, the program uses a matrix library, developed by S. 
White and R. Noack, which links the matrix-matrix and matrix-vector multiplication to corresponding 
BLAS to be more efficient. In order to save memory, only the information for the specific blocks in use 
are loaded, while the rest is stored on disk. For the moment we cannot exceed a number of basis function 
of K = 40 (which uses about 2 Gb of memory) . And right now we are working on rewriting the program 
in order to be able to handle 100 basis functions. 

The program runs under Unix and has mainly be used on workstations such as Linux PC, Dec Alpha, 
HP, SGI and Cray. A tuto rial program for a particle in a box and the matrix library can be found at 
http: / /hedrock. ps.uci.edu/ 
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